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ABSTRACT 

In this paper we investigate the influence of radiative transport on the growth of the 
magnetorotational instability (MRI) in accretion discs. The analysis is performed by 
use of analytical and numerical means. We provide a general dispersion relation to- 
gether with the corresponding eigenfunctions describing the growth rates of small dis- 
turbances on a homogeneous background shear flow. The dispersion relation includes 
compressibility and radiative effects in the flux-limited diffusion approximation. By 
introducing an effective speed of sound, all the effects of radiation transport can be 
subsumed into one single parameter. It can be shown that the growth rates of the ver- 
tical modes - which are the fastest growing ones - are reduced by radiative transport. 
For the case of non-vertical modes, the growth rates may instead be enhanced. We 
quantify the effects of compressibility and radiative diffusion on the growth rates for 
the gas-pressure dominated case. The analytical discussion is supplemented by numer- 
ical simulations, which are also used for a first investigation of the non-linear stage of 
the MRI in gas-pressure dominated accretion discs with radiation transport included. 

Key words: accretion, accretion discs - instabilities - magnetic fields - MHD - 
radiative transfer 



1 INTRODUCTION 

Accretion discs are believed to be present in a variety of as- 
trophysical objects, namely active galactic nuclei, binary X- 
ray sources, cataclysmic variables and protostellar systems. 
The classical problem in the theory of accretion discs is why 
they are accreting, or, stated in other words, what is the 
mechanism which transports angular momentum outwards 
so that matter can spiral inwards? It cannot be ordinary 
molecular viscosity because it is several orders of magni- 
tude too small to account for the observed accretion rates 



(Pringle 1981 Balbus 20031. In order to resolve this dis- 



crepancy, |Shakura fe Syunyaev| ( |1973[ | introduced the phe- 
nomenological a-disc model where a large kinematic viscos- 
ity associated with turbulent motions in the disc is thought 
to drive outward angular momentum transport. What mech- 
anism would make the discs turbulent, however, still remains 
a matter of debate. It is well known that in general shear 
flows are subject to nonlinear hydrodynamical instabilities, 
but it seems that a compelling case for hydrodynamical tur- 
bulence has not been made yet. In fact, the assumption that 
Keplerian discs are hydrodynamically unstable can be chal- 
lenged both on analytical and numerical grounds (see, e.g., 
Hawley et al.|1 999 Balbus & Hawlcy 2006). Since most ac- 
cretion discs consist of ionised gas, it is thus only logical to 
focus on magnetohydrodynamical instabilities instead. 



Nowadays it is the magnetorotational instability (Bal- 
|bus fc HawIey|[T99T| hereafter BH) that appears to be the 
most promising candidate. The instability is local, linear, 
powerful (the growth rates being of the order of the orbital 
period), and it exists even if the initial magnetic field is very 
weak. Both three-dimensional local shearing-box (Hawley 
|et al. 1995 Brandenburg et al.||1995[) a nd global disc sim 



ulations (Armitagc 1998. Hawley 2000 Fromang & Nelson 



2006 1 show that the resulting turbulence indeed transports 
significant amounts of angular momentum outwards, lead- 
ing to values of the a-parameter that are in the range of 
a ~ 10 -3 • • • 10 -2 , depending especially on the initial mag- 
netic field configuration. When compared to observations, it 
seems that numerical simulations generally tend to yield a 
values that are too small, possibly due to numerical effects 
(|King et al.||2007l. The situation now looks rather compli- 



cated, since it has turned out that a does not only depend 
on physical parameters such as viscosity and resistivity, but 
also on details of the numerical setup, namely the numerical 



dissipation and resolution (Pessah et al. 2007 Fromang & 
|Papaloizou|2007 \ . Therefore it is very important to further 
investigate and to quantify the impact of the various physi- 
cal and numerical parameters on the MRI-induced turbulent 
transport. 

Despite the huge amount of numerical and analytical 
work that has been devoted to exploring the MRI in accre- 
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tion discs, comparatively little attention has been paid to 
the issue of radiation transport. Especially when it comes 
to simulations which are global in nature or at least in- 
clude vertical stratification, there are relatively few works 
that include radiation transport (like Turner et. al. 2003, 
Turner 2004 for the radiation pressure-dominated case, Kro- 
lik et. al. 2007, Blaes et. al. 2007 for a disc with comparable 
radiation and gas pressure, and Hirose et. al. 2006 for the 
gas-pressure dominated case). Instead usually an isothermal 
equation of state is used (as, for example, in the works of 
Miller fc Stonel[2000l |Papaloizou fc Nelson||2003l |Fromangl 
fc Nelson|2006 1. Yet from an isothermal model it is not pos- 
sible to derive the true temperature profile, nor can one be 
sure that the resulting vertical structure is correct. If, on 
the other hand, one chooses an adiabatic equation of state, 
but does not include cooling, then the temperature in the 
disc will rise quite rapidly due to turbulent heating, which 
is also not very realistic. Only by including radiative trans- 
fer can we hope to model the energetics of MRI-generated 
turbulence correctly and achieve a quasi-steady turbulent 
state. In this work, we assume that the heat generated by 
the turbulence is transported by radiative diffusion. 

Yet another point is, that incorporating radiative diffu- 
sion might well yield unexpected consequences. For example 
it has turned out recently that the migration of planets in 
a protoplanetary disc is significantly influenced by radiative 
effects (see, e.g., Kley & Crida 2008). It is not unreasonable 
to expect that radiative diffusion might have other impor- 
tant influences, for example when considering the growth 
of planetesimals out of dust in a turbulent disc (this has 
been pointed out recently in |Brandenb urg 2008). For these 
reasons, including radiation transport should be considered 
imperative on the road towards building more realistic mod- 
els of accretion discs. 

The plan of our paper is as follows: In Sec. [2] we de- 
rive the dispersion relation for axisymmetric MRI modes 
in a radiative, gas-pressure dominated accretion disc in the 
flux-limited diffusion approximation and calculate the cor- 
responding eigenfunctions and MRI growth rates. In Sec. [3] 
we describe numerical simulations of single MRI modes us- 
ing two different numerical schemes and demonstrate their 
excellent agreement with the analytical prediction. We also 
investigate the dependence of the saturation level of the MRI 
on the strength of radiative diffusion in a series of 3D local 
shearing box simulations. We conclude in Sec. |4j by dis- 
cussing some consequences of our analysis and providing an 
outlook for further research. In the appendix we relax the 
assumption that gas pressure dominates radiation pressure 
and provide a general dispersion relation, including both the 
radiation dominated and the gas-pressure dominated case. 



As our starting point we take the equations of ideal 
magnetohydrodynamics, namely the continuity equation 



| + v.( P «) = o, 



the equation of motion 

dv „ 1 „ ( B 2 
— +vVv + -V[p+ — 
at p \ 2p 



B VB 

Pop 



and the induction equation 



dB 

~dt 



= V x (v x B). 



(la) 



V<£ = 0, (lb) 



(lc) 



In Eq. \Vo\ , $ denotes the static gravitational potential of 
the central object. It vanishes when linearising the above 
equations, so we do not need to specify it. The other symbols 
have their usual meaning. 

In order to perform a local stability analysis, we con- 
sider a patch of an accretion disc that is small enough so 
that the background density p, pressure p and magnetic field 
B = Btj, (f> + B z z can be taken constant. We do not include a 
radial component in the magnetic field; otherwise the back- 
ground solution would not be time independent because an 
initially radial field would be sheared into a time-dependent 
azimuthal field. For a Keplerian disc, the background flow 
is locally v = — |i?r(/>, with J? the angular orbital fre- 
quency. We impose plane-wave axisymmetric perturbations 
Sp, Sv, SB, Sp proportional to exp[i(fc r r + k z z) + crt], yielding 
the following linearised equations: 



a Sp = — pik ■ Sv, 

1 lfc / 2 . SB B 

a Sv — -S2 Sv r <p H ■ C off Sp H 

2 P \ Po 



(2a) 



\k z B z 
Pop 



SB + 2fl z x Sv = 0, (2b) 



a SB = -lfc ■ SvB + ik z B z Sv 



30 
2a 



S t',. 



(2c) 



where we have introduced the effective sound speed C e s via 
the definition C 2 S = Sp/5p. Note that Eq. ( |2c[ ) ensures that 
the divergence-free constraint k ■ SB = be satisfied. 

After eliminating Sp and SB from the linearised momen- 
tum equation ( 2b l by the use of (2a I and pc|) , the remaining 



system of three equations can be readily solved. By defining 
a wavenumber K via 



K Sv r = k ■ Sv 



(3) 



(which means that K basically constitutes a measure of the 
compressibility of the perturbations), we can write the re- 
sulting dispersion relation in a form that closely resembles 
the incompressible dispersion relation of BH: 



2 LINEAR ANALYSIS 
2.1 Dispersion relation 

We begin our analysis by deriving the MRI dispersion 
relation with radiation transport included in the one- 
temperature flux-limited diffusion approximation (see |Lev^| 
crmore fc Pomraning|1981| |. This means we assume that gas 
and radiation are at the same temperature and that the ra- 
diation energy is small compared to the internal energy of 
the gas, i.e. the disc is gas-pressure dominated. 



k r K 2 



ki 



AQ 2 k 2 v\ z 



2f2av A <t,VA z k z K = 0; (4) 



where k = |fc|, da = B / \/pop denotes the Alfven speed and 
<t 2 = a 2 + k z v% z . The form Q of the dispersion relation ex- 
plicitly demonstrates the dependence on the compressibility. 
The expression for K turns out to be 



K 



a 2 a 2 k r - 2Qakl 



z VA<t>V\ z 



a 2 {a 2 + klCl s ) + o*klv 2 A 



(5) 



A<j> 
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where a is one of the roots of Eq. (J4j . In contrast to what 
one may expect from its definition in Eq. Q, Eq. ^ shows 
K does not couple the modes 5v r and 6v$. In the incom- 
pressible limit C c ff — » oo, we have K — and indeed recover 
the dispersion relation of BH. 

The value of C c ff is determined by the internal energy 
equation. In the one-temperature approximation suitable, 
e.g., for protoplanetary discs, it reads: 

^ + V- (ev) = -pV ■ v - V- F, (6) 

with the radiative flux vector F in the flux-limited diffusion 
approximation 

,VT .., „ ,aT 4 Ac 



F = -eD- 



T 



with D = 4- 



(7) 



Here, k is the opacity, c the speed of light and A the flux lim- 
iter, which depends on the flux-limited diffusion model one 
decides to use. In the linear growth phase, A is a constant, 
which is due to the fact that we consider a uniform back- 
ground. This and the uniformity of the other background 
quantities allows using D as a parameter for our study, which 
will be designated as the "coefficient of radiative diffusion" . 
Linearising Q and using p = (7 — 1) e, as well as the relation 
ST/T = Sp/p — Sp/p we obtain 



aSp = — pC c g ik ■ 5v; 

where C e « is given by the following expression: 

p -y + k 2 D/a 

with 



^eH = 7eff - 



7off 



(8) 



(9) 



p l + k 2 D/a 

We see that in the linear regime, the influence of radiative 
diffusion can be described by introducing an effective adi- 
abatic index 7 c ff. Without radiative diffusion (D — 0), we 
have 7off = 7, and C e a is just the adiabatic gas sound speed. 
With D increasing, the effective adiabatic index 7 e g is re- 
duced towards the isothermal value 7 e g = 1 at D — 00. This 
reduction happens mostly in the regime where the nondi- 
mensional diffusion coefficient k 2 D/a ~ 1, as can be seen 
from Fig. [l] where the effective adiabatic index is plotted as 
a function of D. Thus, in the linear regime, the only effect 
of radiative diffusion is to make the physical situation more 
isothermal. 

In order to aid comparison with the dispersion relations 



derived by other authors, we define (following Blaes & Bal- 
bus|1994||Blaes fc Socrates||2001 \ 

(10a) 
(10b) 
(10c) 



-2 



co - k z (C cS + w A A 



jCj 4 — Q 2 Cu 2 — Af} 2 k 2 v\ z 



and write Eq. Q in an alternative way as 



k z v A(jl v Az (k lo +3k z {2 ) 



k 2 



(lOd) 



= 0. (11) 



In the appropriate limits our dispersion relation is consistent 
with that of other authors, like for example the two fluid 



dispersion relation of Blaes & Balbus (1994) in the limit of 
D — > and zero coupling, or the resistive dispersion relation 
of |Sano fc Miyama| ( |1999[ ) in the limit D — > and zero 
resistvity. 




0.0001 0.01 



100 10000 



Figure 1. Effective adiabatic index f e g as a function of 
fc^ lax D/i7, where fc ma x is the wavenumber of the fastest grow- 
ing mode. In this plot, the azimuthal Alfven speed is taken equal 
to the isothermal sound speed, t)^. = p/p (the global features of 
this plot are not sensitive to this particular choice). 



Blaes & Socrates ( 2001 1 considered the case of a 



radiaton-dominated disc where matter and radiation inter- 
act only by momentum exchange and are, therefore, at dif- 
ferent temperatures, while we consider the opposite limit 
where matter and radiation are closely coupled via emission 
and absorption processes, and are therefore at the same tem- 
perature. In the appendix a more general dispersion relation 
is provided, from which both the Blaes & Socrates dispersion 
relation and the one given in Eq. Q can be derived. 



2.2 Eigenfunctions 

The corresponding eigenfunctions in terms of the radial ve- 
locity perturbation 5v r are given by the following formulae: 



SB 6 =1 



Sp = —ip — • Sv r 
a 

OV z — ; ■ OV r 



Sp = -1/OCefl — ■ 8v r , 

a 

Jn k z B z 

0B r — 1 ■ 0V r , 

a 

k z B z (k 2 b 2 k r K „n 2 \ KB 6 



-3- 



8B Z = — i ^ r ^ z • 8v r 
a 



(12a) 
(12b) 
(12c) 
(12d) 
(12e) 

5v r , 
(12f) 

(12g) 
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(13) 



Using once more the relation ST/T = Sp/p — Sp/p, we can 
derive the following expression for the temperature pertur- 
bation: 

ST _ 7-I 8p 
~T ~ l + k 2 D/a~p' 

As said, radiative diffusion tends to make the situation more 
isothermal, and thus for D — + oo we have ST — > 0, as is to 
be expected. 



2.3 Change of growth rates due to compressibility 
and radiative diffusion 

In order to quantify the effects of a finite compressibility and 
radiative diffusion on the growth rates, we first look at the 
vertical modes (k r — 0) which have the largest growth rates. 
By differentiating the dispersion relation Q with respect to 
C~g, we find for the change of the growth rate a with respect 



to a change in the effective sound speed: 



1 / da 



eff / C. 



! eff=°° 



2k 2 vlJa 
l + 2a 2 /{2 



~2 
l2 V A<j> . 



(14) 



From this result we conclude that for vertical modes the 
effect of a nonzero compressibility leads to a decrease of 
the growth rates in the presence of a nonzero azimuthal 
field. For the fastest growing mode, which, for C e g = oo, 
has a growth rate o" max = 3/4,1? with wavenumber fc max = 
(VT5/4) O/vaz, this means that the change A<r max in the 
growth rate as compared to the incompressible case C c g 
x is approximately: 



^max 



1 V 2 A4> 

"5C7 c 2 ff ' 



(15) 



The relative dampening of the growth rates is, thus, in this 
case of the order of 0(vjL>/Ceff)- 

Curiously, if we consider non-vertical modes (k r 7^ 0) 
and a vanishing azimuthal field, we discover the opposite. 
In this case we find by an analogous analysis: 



da 



dC. 



eff / C, 



~ 2 2/7.2 rft 

a a jk z il 



2k 2 l + 2k 2 a 2 /k 2 n 2 ' 



(16) 



The corresponding shift of the maximum growth rate be- 
comes: 



Aa„ 



27k 2 r 



240fc 2 C 2 S ■ 



(17) 



This effect should not be considered too important, be- 
cause the growth of the MRI will be dominated by the 
fastest growing modes, which are the vertical ones. In both 
cases considered, radiative diffusion contributes a fraction of 
(7P/p)/Ccff — 1 = 7/7cff — 1 to the total shift in the growth 
rates. 

The critical wavelength, which is obtained by setting 
a — in the dispersion relation, is the same as in the non- 
radiative case. This means that radiative diffusion will gen- 
erally not affect the regime in which the MRI may possibly 
operate. One should bear in mind however, that not the lin- 
ear analysis, but only numerical simulations, can provide us 
with the answers if, in a given situation, the turbulence that 
is initiated by the MRI will be sustained for a long period of 
time and how effectively it transports angular momentum. 




(a) k = kz, Ceff »vA 



(c) kr = kz, Ceff « vA 




(b)k = kz,Ceff«vA 



(d)kr = kz,Ceff»vA 



Figure 2. Sketch of the MRI eigenmodes. Plotted is the density 
perturbation 8p (grayscale) and the velocity perturbation Sv (ar- 
rows) in the r-z plane. The plots on the left show the case of a 
vertical mode in the presence of a background azimuthal magnetic 
field; where (a) corresponds to the incompressible case, while (b) 
is for finite compressibility. The plots on the right show the case 
of a nonvertical mode with k r = k z where no azimuthal magnetic 
field is present. Here, (c) corresponds to the case of maximum 
compressibility and (d) to the incompressible case. 



For example, recent results show that in zero-net flux shear- 
ing box calculations this depends critically on the value of 
the magnetic Prandtl number ( Fromang et al.|2007 l. 

To summarise, radiative diffusion changes nothing fun- 
damental regarding the growth of the MRI. This is because 
the radiative effects enter only via the effective adiabatic 
index y e g and, thus, change the compressibility, but noth- 
ing else. This statement, that the growth of the MRI under 
the influence of radiative diffusion can be characterised by 
an effective sound speed, remains true also if we drop the 
one-temperature approximation (cf. the appendix). 

In gas-pressure dominated discs, the growth rates will 
never be dramatically changed by radiative diffusion, since 
the variation in C 2 S covers only a factor of 7. This is differ- 
ent from the situation encountered in a radiation dominated 
disc, where the analog of Cfg may vary greatly. In the pres- 
ence of a strong enough azimuthal field, the growth rates 
may thus be severely reduced (cf. the appendix and the pa- 
per of |Turner et al.|2002| . 

In order to understand the change of the growth rates 
due to radiative diffusion in a qualitative manner, let us con- 
sider two cases: First, the case of a vertical mode (fe = k z e z ) 
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in the presence of an azimuthal field. In the incompress- 
ible limit (C e ff — > oo), the motion of the fluid is confined 
to the plane perpendicular to the perturbation wavevector, 
so that 8v z — (Fig. [2] (a)). When the compressibility is 
nonzero, the Lorentz force due to the azimuthal magnetic 
field causes a fluid flow in the vertical direction (Fig.[2](b)). 
The higher the compressibility (and the smaller therefore 
the gas pressure), the more vertical will the velocity vector 
become. This in turn makes the buildup of the magnetic 
field less effective, resulting in a smaller growth rate. Ra- 
diative diffusion increases the compressibility, and therefore 
decreases the growth rate. 

Next consider the case of a nonvertical mode (k r 7^ 0), 
with no background azimuthal field. If we now start with 
the limit of maximum compressibility (C e a = 0), we have 
Sv z — (Fig. [2] (c)), since the z-component of the Lorentz 
force vanishes. When now increasing C e g, the gas pressure 
will act to push the velocity vector into the plane perpen- 
dicular to the perturbation wavevector, this effect becom- 
ing stronger and stronger as the compressibility is further 
decreased (Fig. [5] (d)). Therefore, this time increasing the 
compressibility makes the buildup of the magnetic field more 
effective, which means that the growth rate increases due to 
radiative diffusion. 

When considering the general case of a nonvertical 
mode in the presence of a nonzero azimuthal field, both ef- 
fects are present and the result will be either an increase 
or a decrease of the growth rate, depending on the strength 
of the azimuthal field and the direction of the perturbation 
wavevector. 



3 NUMERICAL SIMULATIONS 

In this section we investigate the growth of the MRI using 
the appropriate numerical tools. Here we especially aim at 
a verification of these numerical tools by a comparison to 
the analytical results derived above. After discussing the 
numerical scheme we will present the results of the numerical 
simulations of the MRI. 



3.1 Numerical approach 

The investigation of the analytical results using a numerical 
method is clearly not very demanding in this case owing to 
the fact that the analytical model is restricted to the lin- 
ear growth phase. Here we decided, however, for a scheme, 
which is well adapted to high Mach-number turbulence sim- 
ulations, because we are also interested in more general ac- 
cretion disc simulations including radiation transport. One 
of the most important constraints for a numerical scheme 
used for compressible turbulence simulations is the correct 
handling of sharp discontinuities. Here we use a second or- 
der finite volume scheme based on the work by |Kurganov] 
et al. (2001 I for the hyperbolic part of the system of equa- 



e.g., |Balsara fc Spicer|p99| |Londrillo fc Del Zanna||2000 1 



tions. This is a central conservative scheme for the solution 
of equations of the form: 



dt 



+ V ■ F(m) = 0. 



(18) 



It does not require a Riemann solver and a characteris- 
tic decomposition. This scheme was extended by a con- 
strained transport description for the magnetic field (see, 



by which we assure the solenoidality of the magnetic field. 
This method uses the hyperbolic fluxes to compute electric 
fields on a staggered grid. These are then used to evolve 
the magnetic induction, the components of which are also 
given on a (different) staggered grid. The stability of the 
code and its capability to resolve steep gradients without 
introducing artificial oscillations have been proven, e.g., in 



Kissmann| fl2006| ). 

Being interested in the evolution of the gas under the 
influence of radiation transport we also have to deal with the 
corresponding source terms to the energy equation. These 
were implemented using two alternative approaches. In the 
first approach we included the diffusive radiation transport 
explicitly in the scheme, whereas we used an implicit solver 
for the second method. This dual approach is motivated by 
the fact that we are interested in being able to cross-check 
the two different methods. 

For the implicit implementation of the source terms we 
split the evolution of the radiation from the evolution of the 
hyperbolic part of the system of equations: 



+ 



(*) 1 



hyp 



(19) 



(20) 



where we evolve the radiation transport step using a parallel 
oj-Jacobi matrix solver. The oj-Jacobi method is very similar 
to the method of successive overrelaxation (SOR). It has the 
advantage that it is simple to implement and very easy to 
parallelise. We have also a multigrid method available, in 
case that the w-Jacobi solver should prove insufficient for 
high-resolution simulations. 

Having an implicit solver for the radiation transport 
is of special importance with regard to future applications. 
That is for high resolution simulations with strong radiation 
transport the time-step of an explicit scheme would become 
prohibitively small, whereas an implicit solver would only be 
limited by the Courant condition (see Courant et al.|[l928 \ 
for the hyperbolic part of the system. 

Our second implementation is nonetheless of explicit 
form, because we wanted to have a completely different 
method at hand to be able to compare it to the implicit 
scheme for more general simulations for which no analyti- 
cal solution is available. For this purpose we extended the 
flux-function for the energy equation by: 



p Rad — _ £ n 
T 



dT dT dT 
x ay 



dz 



(21) 



Since the flux functions are only needed at the cell faces the 
occurring gradients are easily computed as can be seen from 
the example of the ^-component: 



riRad 

- x,i+l/2,j,k 



T 



T, 



Ax 



This scheme proved to be stable for 
At < 



Ax 2 pT/e 



(22) 



(23) 



8D 7- 1 

In the following section we will show that both implemen- 
tations for the radiation transport are well-behaved, yield- 
ing correct result when applied to the analytical problem at 
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Figure 3. MRI growth rates obtained with the implicit solver. 
Plotted are the analytical solutions for various cases, together 
with corresponding data points obtained by numerical simula- 
tions. The curves in the main part of the graphic are for vertical 
modes (k r = 0), where the growth rates are reduced by radiative 
diffusion. From top to bottom, they represent the following cases, 
= D = 0; = 20B Z ,D = 1(T 5 ; B^ = 20B Z ,D = 10" 4 ; 
B^ = 20 B Z ,D = 1. The vertical magnetic field in terms of the 
plasma beta, fi = 2pop/B 2 , is chosen such that /3 = 400, which 
means that VAz/yP/p — 0.1. The inset shows simulation results 
for nonvertical modes with no azimuthal magnetic field, where 
radiative diffusion increases the growth rates. We chose k r = k z , 
B^ = and B z such that Vj\ z = \/p/ p. From bottom to top, 
we have D = 0, D = 10 -5 , D = 10 -4 , D = 1. Although the 
difference in the growth rates is only small for the latter case, it 
is accurately reproduced by our numerical code. 



hand. Thus, they can be used for future simulations of more 
complex radiation transport phenomena. 



3.2 Growth rates 

For our numerical simulations, we use the shearing box ap- 
proximation (see, e.g., Hawley et al.|1995 |. In order to com- 
pare the analytical growth rates of the MRI to the results 
from our numerical code, we perform simulations of single 
MRI eigenmodes. Since only variations in the radial and in 
the z-directions have to be taken into account, it is sufficient 
to use a 2D model. We set up the problem using the eigen- 
functions as given by (12 1. The size of the computational 
domain is [-0.5,0.5] x [-0.5,0.5] (x- and z-direction respec- 
tively) with a typical resolution of 64x256 cells. Initially we 
set fl = 10~ 3 , density p — 1 and pressure p — 10~ 6 . The 
simulations were performed for different wavenumbers and 
for various strengths of magnetic field and radiative diffusion 
coefficient (for the explicit values see in the figures) . 

The results are plotted in Figs. [3] and [4] for the im- 
plicit and the explicit implementation of the source terms, 
respectively, together with the corresponding analytical so- 
lution. In the simulations with the explicit scheme a con- 
stant diffusion coefficient was used [cf. Eq. (21 1], while in 




k z v Az /tl 

Figure 4. Comparison of the numerical growth rates (black 
crosses) to the analytical model results (full lines). Here we show 
the data from the simulations with the explicit numerical im- 
plementation of the radiation transport terms. All simulations 
have B^ = 25B Z and a vertical magnetic field corresponding to a 
plasma beta j3 = 400. The curve for D = oo has been inserted to 
show the maximum possible change in the growth rates. 



the simulations with the implicit scheme, the full radiation 
transport was done as prescribed by Eq. |7|), with A = 1/3 
and constant opacity k. The data points nicely match the 
analytical prediction, thereby confirming both the validity 
of the numerical scheme and the analytical calculation. In 
particular it is evident that both implementations of the 
radiation transport are in good agreement to the analytical 
predictions. Therefore, we can use both methods to simulate 
more complex situations in the future giving us the good op- 
portunity to have two different methods to be applied to the 
same problem. Thus, it will be easier to decide if observed 
effects might be due to some numerical problem. 



3.3 Saturation level 

In order to study the impact of radiation transport on the 
saturation level of the MRI, we switch to a 3D representation 
of the shearing box with a computational domain [-0.5,0.5] 
x [0,4] x [-0.5,0.5]. Concerning the values of p,p, Q we use a 
similar setup as in the 2D simulations described in the pre- 
vious section. We choose an initially vertical magnetic field 
corresponding to a plasma beta of j3 = 400 in all simula- 
tions. Instead of prescribing eigenfunctions, we initialise the 
problem with random velocity and pressure perturbations 
of order 10" 6 . 

The strength of the angular momentum transport is 
described by the a-parameter ( Shakura fc Syunyaev|[l973 1 . 
In our simulations, it is measured according to the following 
prescription: 



q= (T R + T M )/po, 



(24) 



where po is the initial gas pressure, Tr = p8v r Sv,p denotes 
the Reynolds stress and Tm = —SBrSB,/, is the Maxwell 
stress. The angular brackets (■ • ■ ) indicate a spatial average. 
We also use the dimensionless stresses qr and om defined 
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Table 1. Time- and space-averaged quantities from 3D shearing 
box simulations performed with our code. In all of the simulations 
B,p = 0. In the first three lines, only the resolution is varied, 
while the last three lines show the dependence on the strength of 
radiative diffusion. Averages were taken over 10 to 40 orbits. 



Resolution 


D 


-Emag/P0 


a 


a M 


OR 


32x64x32 


1 


0.400 


0.278 


0.248 


0.044 


128x256x128 


1 


0.505 


0.299 


0.279 


0.036 


64x128x64 


1 


0.337 


0.226 


0.195 


0.031 


64x128x64 





1.007 


0.719 


0.599 


0.119 


64x128x64 


10~ 5 


0.358 


0.244 


0.210 


0.034 


64x128x64 


1 


0.337 


0.226 


0.195 


0.031 



as 

or = <Tr)/Po, a M = (T M )/p - (25) 

Note that with this definitions in general we will have 
a 7^ or + om- Due to the fact that in our simulations the 
magnetic field posesses a net flux, the resulting a is quite 
large, a > 0.1. 

In order to check if there exists a trend for the turbu- 
lent activity with different resolution, we performed three 
runs where the resolution was successively doubled. Fig. [5] 
shows the time evolution of the magnetic energy and the 
a-parameter for this runs. In Table [I] long-term averages 
of this quantities as well as of the Maxwell and Reynold 
stresses can be found. The results displayed in Fig. [5] and 
Tab. [I] indicate no net trend for the change of the saturation 
level with resolution (the first three lines in the table). 

As we have seen, the analytical calculation predicts that 
radiative diffusion modifies the MRI growth rates. The ques- 
tion if the saturation level of the MRI-induced turbulence 
is also changed by radiative effects can only be answered by 
numerical simulations. We have performed simulations with 
varying radiative diffusion coefficient, where we picked a res- 
olution of 64x128x64. The results show a moderate decrease 
of the turbulent activity with inreasing radiative diffusion 
coefficient D, see last 3 entries in Tab. [T] A qualitative ex- 
planation for this is given below. 

We have also performed additional simulations with the 
ZEUS code in order to check this result using a different nu- 
merical scheme. We used the same setup concerning p,p,0 as 
in the simulations done with our code. The strength of the 
constant vertical magnetic field was taken as corresponding 
to a plasma beta of /3 = 800. We performed simulations 
with and without azimuthal magnetic field using dif- 
ferent resolutions (see Table [2J. In Fig. [6] the magnetic en- 
ergy is plotted for the high resolution simulations. As has 
already been stated in the introduction, it is known that 
the outcome of an MRI simulation does depend on numeri- 
cal issues such as the resolution or the numerical magnetic 
Prandtl number of the code ( |Fromang fc Papaloizou|2007[ ). 
Indeed, if we compare for example the simulations in lines 
3 and 4 of Tab. [l] with the corresponding simulations done 
with the ZEUS code, lines 5 and 6 of Tab. [2] the ratio of 
the saturation levels in the case of the ZEUS code turns 
out to be about a factor of two smaller than in the simu- 
lations done with our code. We have checked that this fact 
remains true even if we choose the same initial plasma beta 



Table 2. 3D shearing-box simulations performed with ZEUS. In 
some of the simulations, a non-zero constant was used, leading 
to a higher saturation level. Averages were taken from 10 to 40 
orbits. 



Resolution 


D 


B z 


-Emag/P0 


a 


Tyi 


7r 


32x64x32 








0.398 


0.248 


0.210 


0.038 


32x64x32 


1 





0.245 


0.146 


0.123 


0.023 


32x64x32 





20 


2.557 


1.125 


0.983 


0.141 


32x64x32 


1 


20 


1.542 


0.577 


0.498 


0.079 


64x128x64 








0.552 


0.335 


0.288 


0.047 


64x128x64 


1 





0.392 


0.218 


0.187 


0.031 


64x128x64 





20 


2.474 


1.137 


1.006 


0.131 


64x128x64 


1 


20 


2.315 


1.075 


0.938 


0.137 



as in the runs with our code. Despite these issues, all runs 
unambigously show that the saturation level decreases with 
increasing radiative diffusion. 

A plausible explanation for the phenomenon of the re- 
duction of the saturation level due to radiative diffusion can 
be given as follows: In the case of an initial magnetic field 
with nonzero net flux, the turbulence in the saturated state 
is still partially pumped by the two-channel mode (the verti- 
cal mode with the longest wavelength that fits into the box). 
The channel mode reappears every few orbits; a manifesta- 
tion of this are the oscillations that can be seen in the time 
evolution of the magnetic energy and the cv-parameter (see, 
for example [Sano fc Inutsuka|200l| . In Sec.|2~3| we have seen 
that radiative transport tends to decrease the growth rates 
of vertical modes by effectively increasing the compressibil- 
ity of the fluid and thus making it easier for the Lorentz force 
to push the velocity vectors out of the r-<f> plane. Therefore it 
is to be expected that radiative transport, by decreasing the 
growth rate of the recurrent channel mode, acts to reduce 
the turbulent activity. 

3.4 Temperature distribution 

For the full 3D simulation the influence of the radiation 
transport can also be visualised by the temperature distribu- 
tion. This is shown in Fig.[7]for a similar setup as described 
in the previous section. Here we used an extent of the com- 
putational domain of [-0.5,0.5] x [0,6] x [-0.5,0.5] with the 
same spatial resolution of 64x128x64. The ^-component of 
the magnetic induction was initialised with /3 = 800, and 
we chose B^, — 20B Z . We followed the evolution for several 
orbits. In Fig. [8] we show the distribution in the density - 
thermal energy plane for one simulation without and one 
with radiative transport. Apparently the plasma gets more 
isothermal when radiation transport is present - for a fully 
isothermal simulation the thermal energy would depend lin- 
early on the density, thus, yielding a line in this plot. This 
fact also becomes apparent, when comparing the results for 
the temperature distribution function. 

This is done in Fig. [7] for the same simulations. Here, 
the temperatures are normalised to the initial temperature 
(that is, initially we had a delta-peak at T = 1). The shift 
of the distribution functions with regard to the initial tem- 
perature results from the heating by (numerical) viscosity 
and resistivity. Obviously the distribution for the disc with 
radiation transport is considerably narrower. Note, however, 



10 c 



10' 



i<r- i 



10-* 




32x64x32 

64x128x64 

128x256x128 



32x64x32 

64x128x64 

128x256x128 



10 



20 
t [orbits] 



30 



40 




20 
t [orbits] 



Figure 5. Resolution dependence of the turbulent activity on the numerical resolution. The left and right panels show the time evolution 
of the perturbed magnetic energy and the alpha parameter, respectively. In case of the alpha parameter the curves have been smoothed 
over 1 orbit. When considering long-term averages of these quantities, there appears to be no trend with varying resolution. 





Figure 7. Temperature distribution for two 3D simulations without radiation transport (left) and with radiative transport, D = 4- 10 -5 
(right). In both plots, the distribution to the left corresponds to t = 3.5 orbits, while the right corresponds to t = 4.1 orbits. In the case 
with radiative transport, the width of the distribution function is smaller and the heating occurs less rapidly. 



that also for this case the temperature after several orbits 
differs markedly from the initial temperature. This conlu- 
sion is also strengthend by the fact that both distributions 
are centered at nearly the same temperature. This is due to 
the fact that the gas is still heated by viscosity and resistiv- 
ity. In contrast to the case without radiation transport the 
local temperature enhancements are smoothed out by the 
radiation transport. 



4 DISCUSSION 

We analysed the influence of radiation transport on the mag- 
netorotational instability. We started by investigating the 
linear growth of small disturbances on an underlying con- 
stant shear flow. The system under investigation is a mag- 
netised shearing box, where we explicitly retained the com- 
pressibility of the gas. Additionally we took radiation trans- 
port via flux-limited diffusion into account. For this system 
we derived a general dispersion relation, which in the appro- 
priate limits can be reduced to several dispersion relations 
known from the literature. 
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1e-12 



Figure 6. 3D shearing box simulations performed with ZEUS 
at resolution 64x128x64. Shown is the perturbed magnetic en- 
ergy normalised to the initial gas pressure. The solid, dashed, 
dot-dashed and dotted curves correspond to the combinations 



0/D 



o, .a* 



0/D 



1, B A 



20B Z /D = and 



Bff, = 20B Z /D = 1, respectively. In the linear phase, the obtained 
growth rates are only a few percent smaller than the growth rate 
of the fastest growing mode that fits in the box, which means that 
the initial growth is dominated by the fastest growing mode, as 
is to be expected. 



Here, we placed special emphasis on the one- 
temperature flux-limited diffusion approximation suited for 
accretion discs not dominated by the radiation. The corre- 
sponding dispersion relation can also be derived from the 
general relation given in the appendix. From this we com- 
puted the growth rates for the linear case, which we then 
also compared to the results from numerical simulations. 
The excellent agreement between simulations and analytical 



prediction not only confirmed the accuracy of the employed 
numerical schemes but also the correctness of the analytical 
solution. 

This analysis was then further extended by numerical 
simulations, which also took the non-linear phase of the in- 
stability into account. These simulations were done for a 
more general case, where several modes of the instability 
where excited in a fully three dimensional setup. There we 
concentrated on the resulting saturation levels of the tur- 
bulence and on the temperature distribution, which is also 
influenced by the radiation transport. 

From our analysis it is apparent that the influence of the 
radiation transport is a rich phenomenon. It not only weak- 
ens the magnetorotational instability, but can also enhance 
it in a certain parameter regime. This shows that radiation 
transport - even in the one-temperature flux-limited diffu- 
sion limit - does not simply act as a dissipative process like 
viscosity and resistivity. Another difference between these 
processes is that the latter can completely suppress the in- 
stability, whereas radiation transport can only weaken it. 
This is due to the fact that a diffusive radiation transport 
only smoothes the temperature, whereas the other dissipa- 
tive processes directly act onto either velocity or magnetic 
field fluctuations. Thus, we have, for strong radiation trans- 
port, approximately an isothermal state of the plasma. This, 
however, still differs from the use of an isothermal equation 
of state with the difference arising from the fact that when 
using an energy equation with 7 differing from one, we still 
have heating due to dissipative processes (either physical or 
numerical). Thus, the effect of radiation transport is just to 
smooth out the local peaks of heat appearing due to local 
dissipation. 

This picture will probably become different, when in- 
vestigating the MRI in an open box instead of the periodic 
shearing box. The future task will, thus, be the investigation 
of a stratified shearing box, where we can allow for outgoing 
radiation in the halo of the disc. This way even a nearly 
isothermal core of the disc might occur, which might then 
be compared to simulations using an isothermal equation of 
state. Here, however, we have restricted ourselves to the in- 
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vestigation of the local shearing box, since only this can be 
approached by analytical means. 
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APPENDIX A: GENERAL DISPERSION 
RELATION 

In this appendix, we provide the general dispersion relation 
by including the radiation energy equation and dropping the 
one-temperature approximation. This means that in the mo- 



mentum equation ( lb l, we have to replace Vp by Vp + \X7E, 
where E denotes the radiation energy density. As in Eq. |7| 
A denotes the flux limiter, which is constant in the linear 
growth phase due to the constant background. The effect of 
this is that in the linearised momentum equation (2b I, we 
should now define the effective sound speed C e f{ as 



2 _ Sp SE 

'-'off = j r 

Sp op 



(Al) 



To determine it, we need the internal energy equation and 
the radiation energy equation: 



dE 

~dt 



+ V ■ (Ev) 



-pV • v - K P p(47tB - cE), (A2a) 
E 

-— V • v + n P p(4nB - cE) - V • F. 

(A2b) 



Here, B denotes the Planck function and ftp the Planck 
mean opacity. The terms containing ftp describe the coupling 
of the gas to the radiation via emission and absorption. The 
radiative flux vector now reads 

Ac 



F = — —X7E, 
kfP 

where kf is the flux-mean total opacity. We define 



Qp 



ftppc 



and qf 



k 2 Ac 



(A3) 



(A4) 



The dimensionless quantity ap gives the rate at which ra- 
diation and matter are thermally coupled with respect to 
the growth rate of the instability while ap constitutes a 
non-dimensional diffusion coefficient. We assume that in the 
equilibrium the gas temperature T s and radiation temper- 
ature T r defined by E = aT 4 are equal: T g — T T = T and 
4 ^p- = aT 4 . Now the linearised equations corresponding to 
Eqs. iA2\ become: 



Sp = 7 — Sp — (7 
P 



l)a P E I I 
ST a 



AE 

8E = —Sp + apE ( 4 

3p V T 



8T± 
T 
SE 

1e 



SE 

~E 

- ap 



SE; 



(A5) 
(A6) 



From these linearised equations we find after some algebra 

Sp _ (1 + a P + a F ) (7 + |f ) - |f (1 - 3a P )(l + a F ) p 
Sp' 



1 + ap + «f + 4ap — (1 + a F ) 



5_E 
Jp~ 



1 + 4ap (7 — 1 + 



E 



1 + ap + a F + 4ap - 



a F ) p 



P 

(A7a) 
(A7b) 



In conclusion, let us look at the limits of gas-pressure domi- 
nated and radiation-pressure dominated discs, respectively. 
In the limit ap — + 00, E/e — > 0, suitable for a gas-pressure 
dominated disc, we get 



Sp 

T P 



7 + 4-a F p 
l + 4fa F p' 



SE 
Jp 



0; 



(A8) 



and in this way recover our Eqs. Q. For the case of a radi- 
ation dominated disc where ap — > 0, 

cls = 1 j> + AXEjZp 
p 1 + a F 

When using this value for C e ff, and setting A = 1/3, Eq. (Ill 



becomes identical to the Blaes & Socrates (2001 1 dispersion 



relation. The general dispersion relation provided in this ap- 
pendix, thus, encompasses both the case of a gas-pressure 
dominated disc, such as a protoplanetary disc, and the oppo- 
site extreme of a radiation-pressure dominated system like 
an accretion disc around a black hole. 
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